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O l' We address the issue of how to identify the equations of a largely unknown chaotic system from 

' knowledge of its state evolution. The technique can be applied to the estimation of parameters that 

drift slowly with time. To accomplish this, we propose an adaptive strategy that aims at synchro- 
nizing the unknown real system with another system whose parameters are adaptively evolved to 
' converge on those of the real one. Our proposed strategy is tested to identify the equations of the 

Lorenz, and the Rossler systems. We also consider the effects of measurement noise and of deviation 
, — of our fitting model from consistency with the true dynamics. 

B 

' O ' In this paper, we present a new technique to estimate the equations of an unknown chaotic system. 
^ We assume that suitable but rather minimal information is available about the form of the system 
equations and we use the observed dynamics to estimate an approximation to the unknown system. We 
. envision that our findings could be useful in situations in which the equations governing the dynamics 

(N 

of a given unknown chaotic system are to be identified or in situations where system parameters drift 



O 
O 



with time. 

I. INTRODUCTION 

Systems of nonlinear differential equations are often used to study the dynamics of real world systems. For instance, 
the Hindmarsh-Rose equation ]| is widely considered to be a reasonable model of the firing/bursting behavior in real 
neurons, and has been shown to be chaotic for a certain range of its parameters. Another example is that of the 
Chua system 2|, |3| that models the dynamics of a simple electronic circuit for which the emergence of chaos has been 
observed both in simulations and in experiments. A general question that often arises is whether such a system of 
differential equations replicates the dynamics of a given real system that it is meant to model. 

In a recent paper [4] , it has been shown that it is possible to synchronize a real experimental system with a system 



of differential equations simulated on a computer by coupling the two. This is only possible if the system simulated 
on the computer does a sufficiently good job of quantitatively replicating the dynamics of the experimental system 
5|. An interesting implication of this is that it is possible to quantitatively evaluate the degree to which a given 
mathematical model replicates the dynamics of a real system via its capability to synchronize with the real system 
when the two are coupled. Here we adopt a different but related application of synchronism of coupled chaotic systems. 
In particular, we consider that the real system is largely unknown, and we propose an adaptive strategy that, by using 
synchronization, is able to obtain a set of (nonlinear) differential equations that describes it. Furthermore, in the case 
that the fitting function basis for our model is not consistent with the true system dynamics, we will try to answer 
the question of how well the obtained model is able to forecast the true system future behavior. 

Refs. P, Ol have outlined the connection between the problem of synchronization of dynamical systems and the 
problem of the design of an observer to reconstruct the state of an unknown system. The idea of using synchronization 
or control for parameter and model identification has previously been presented in P] and subsequently in {q. 
Some recent papers have appeared on reconstructing the state and parameters of the Chua system 
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Here we 



will address the issue of how to identify the equations of a largely unknown generic system from knowledge of its state 
evolution. In Sec. II, the problem is introduced and an adaptive strategy, based on the minimization of appropriately 
defined potentials (see e.g., [M]), is presented. In Sec. Ill, our proposed strategy is tested to identify the equations 
of the Rossler, and the Lorenz systems. In Sec. IV, we will take into account the effects of measurement noise. In 
Sec. V, we will address the case that the fitting function basis for our model is not consistent with the true system 
equations and for this case, we will study the capability of the obtained approximate model to forecast the evolution 
of the true system. 



II. FORMULATION 

We consider the general example, x ~ F{x.), with x — {xi,X2, ...,x„)-^ and F(x) — [/i (x) , /2 (x) , /„(x)]'^, where 

n n n 

fi{-x) = ^ ^ ai-jkXjXk + ^ hjx-j + Ci. (1) 
j=i k=j j=i 

That is, /i(x) is a degree 2 polynomial. Note that the specification of i^(x) from involves M ~ [{n'^ +n) /2+n + l]n 
parameters. For example, in the case of the Rossler system (to be used later in our numerical experiments), n = 3, 
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and 



F(x) = 



-X2 - X3 

xi + 0.165a;2 
0.2 + (xi - 10)a;3 



(2) 



and all the coefficients {aijk,bij, and c^) are zero, except for 0313 = 1, 612 — —1, 613 = —1, 621 — 1, &22 = 0.165, 
&33 = —10, C3 = 0.2. The Lorenz system could be recast in similar terms. The Hindmarsh-Rose system, has an 
analogous structure, but with higher order powers of the state variables as well (one term is third order). 

In this paper, we will be mainly concerned with the case that the chosen form for our model system is a degree 2 
polynomial, that is as in Eq. [191 Nonetheless the extension to the case of higher order polynomials or to a different 
fitting functions basis is straightforward. 

We assume that the exact system (2) is unknown, but it is known that the system is of the form x = -F(x), with x 
n-dimensional, and -F(x) expressible, or approximately expressible, in terms of degree 2 polynomials in x as discussed 
above. It is then appropriate to try to model the dynamics of the true system by x' = F'(x'), with 



x'i — fi{x[,X2, ■■■,x'n) — ^ ^ a'^ji^x'jX'i^ + ^ KjX'j + c'i- 
j=i k=j j=i 



(3) 



Our goal is then to obtain good estimates a'^j^., b'^j, and of the true coefficients aijk, bij, and c^. We assume 
that, although -F(x) is unknown, we do have access to good measurements of the evolving 'experimental' system. To 
accomplish our goal, we envision coupling the true system to the model system (3). Our approach will be to formulate 
an adaptive procedure which adjusts the model coefficients a^jj,, fo'y, and c[ in such a way as to achieve synchrony. 
It is then hoped that, when this has been accomplished, the model coefficients will be a good approximation to the 
corresponding coefficients of the real system. 

We perform a one way diffusive coupling from the true system to the model, as follows 



x' =F'(x') + r(i/(x) -i?(x')). 



(4) 



The quantity H is in general an m < n vector of m observable scalar quantities that are assumed to be known 
functions of the system state x(f). F is an m x n constant coupling matrix. For our numerical experiments we will 
assume that: (i) i?(x) — x; (ii) F — 7/„, where 7 is a scalar quantity, and /„ is the identity matrix of dimension n. 



Note that our strategy requires that, when the system parameters are correctly identified, i.e., a'^-j, = aijk, b'^j = bij, 
c[ = Ci, ^ belongs to the range for which the synchronized solution x(i) — x'(i) is stable with respect to infinitesimal 



perturbations 



15|. 



We now introduce the following potentials, 

=< [i^ >^, i = l,2,...,7z, (5) 

where < G{t) >^ denotes the sliding exponential average v f e-'^(t-t')G{t')dt'. Note that is a function 
of time and also of the coefficients a'^j^., b[p and c[. Also note that if j:2(i), ...,a;„(i) are chaotic, then 

/i(xi(i),a::2(i),...,a;„(i)),/2(a;i(t),a;2(t),...,a;„(t)),...,/„(a::i(i),X2(i),---,2:„(i)) vary chaotically, as weU. Furthermore, 
we point out that the exponential averaging operation < G(t) >i, is the same as low-pass filtering of G(t), using a 
first order filter of bandwidth v. The potential (5) satisfies ^'i > and has a minimum value of zero. A sufficient 
condition for the potentials in ([5]) to be zero is 

Xi{t)=x\{t), i = l,2,...,n, (6) 

and 

'^'ijk — C'ijk, i, j = 1, ••■7 n;k = j, ...n, 

Kj^bi-i, i,j = l,...,n, (7) 

C^, 1 1,..,71. 

Equation ^ corresponds to synchronization. Equation ([7]) corresponds to a correct identification of the (unknown) 
system parameters. Due to the chaotic nature of the true system, it should typically be the case that this sufficient 
condition for the potential to be zero is also necessary. 

We note that in the case in which our proposed model form is consistent with the dynamics of the true system 
(e.g., they are both expressed by a degree two-polynomial as assumed in Eqs. p9p and there is only one possible 
choice of a'ijk,b[^, and that minimizes the potentials ([5]), namely , a'-^j. = aijk, b^^ — bij, c[ = This follows from 
the principle that two polynomials are equal only if all their coefficients are equal and from the fact that ([5]) implies 
a one to one correspondence between the model state variables and those of the true system. 
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We propose to adaptively evolve the estimates of the parameters a'ijk^Kjy^i time, according to the following 
gradient descent relations: 

%^=-/3.^, (8a) 
at oa'ijk 

-l^-'^'db^' ^^^^ 

dc\{t) d^i 

0a,0b,/3c > 0. Our hope is that a-^j,(t), b[j{t), and c^(f) will converge under this evolution to the true parameter 

values, atjk,bij,Ci. 

First we consider (8a). Let {fi)jk denote f-{x[,x'2, ■■■,x'^) evaluated at a'^j^ = 0, 



fl{x[,X2, ...,x'J = a'ijkx'jx'k + {f'i)jk, (9) 



Substituting this into the right hand side of Eq. (8a), we obtain, 



= -2/3„ < a[,,x',\',' + {f^ikx'A - x,x',x'^ >. . (10) 
Similarly letting (/j')j denote f-{x[,x'2..., ,x'^) evaluated at b'^j = 0, Eq. (8b) gives 

= -2/3, < bijx'/ + (/;),x;. - x,x', >. . (11) 

Finally, we consider relation (8c) with (//) denoting fi{x[,X2, ■■■,x'^) evaluated at c- = 0. Then 

= -2/3, <4 + (/;) -Xi>^. (12) 



dt 



In this paper we consider the case where Pa,b,c are very large. For this situation the solutions a'^j^, b'^j, c'^ rapidly 
converge to the minimum of the potentials (which is zero). Thus we can set the the averages < ... >i, on the right 
hand side of Eqs. (10)-(12) to zero. We further consider that the average < ... >y is performed over a time scale 
that is much larger than the time scale Tg for variation of x(t), in which case a'ijkit), b^j{t), c^{t) vary slowly 
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compared to x(t). Under these conditions, (8), (10), (11), and (12) then yield 



n n n 



j = l,...,n; k 



j,...,n, 




(13a) 



n n n 



'Jiim, < ^i^;^^;^- >y + 6^; < x'lx'j >^ < x'j >^— < Xix'j > 



J 




(13b) 



n n n 



(13c) 




Equations (13) constitute a system of M = [(n^/2) + (3n/2) + l]n Unear equations for the M quantities a'-j/^, b^j, 
and C-. The coefficients of the quantities to be solved for, as well as the driving factors on the right hand sides of 
Eqs. (13), are all of the form of an average < (...) where for the coefficients of the unknown (...) is a product of 
x' terms, while for the driving factors (...) involve the time derivative x of the observed experimental system state. 
In practice, it is inconvenient to explicitly calculate the integrals for these quantities, in terms of the form. 



and obtain I{t) as a function of time by solving (15). Thus our adaptive system for finding estimates of the quantities, 
dijkjbijjCi, is Eq. (4) for x'(t), Eqs. (13) for a'-ji.{t), b[j{t), and c[{t), where the various terms in (13) that are of the 
form I{t) =< G{t) >v obtained by integrating (15). If our procedure works, the time evolutions of a'^jhi^), b[j{t), c[{t) 
will converge to aijk, hj, Ci with increasing t. 

In this paper, we focus on a case in which the true system equations are linear in the unknown coefficients (e.g., 
they are expressible in forms of given degree polynomials as in (1)). Under this assumption, the minimization of the 
potentials can be achieved by solving a system of linear equations as in (13), which in practice, is simpler than solving 
the system of differential equations (10-12). We note, however, that, our strategy can also be employed, when the 
true system equations (1) are nonlinear in the unknown coefficients. For such a case, it may be impossible to obtain 




(14) 



at every time step. Instead we shall use the fact that I{t) satisfies the differential equation. 



(15) 
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a simple unique solution for the unknown coefficients as in (13), and integration of the gradient descent differential 
equations (e.g., Eqs. (10-12)) is a potentially useful approach. 

III. NUMERICAL EXPERIMENTS 
A. An experiment with the Rossler system 

We now present numerical experiments testing the above strategy for the example in which the unknown system 
is the Rossler system described by Eq. ([2]). The system in ([2]) is evolved starting with a random initial condition on 
the attractor x{0),y{0), z(0), while the system in ^ is evolved starting from a perturbed initial condition, 

x[{0) = xi{0) + pie,, 4(0) =X2(0)+p2ey, 4(0) = 2:3(0) + psie.l, (16) 

where ex,ey and ez are zero- mean independent random numbers of unit variance drawn from a normal distribution; 
Pi — 7.45, p2 — 7.08, p3 — 4.25 are the standard deviations of the time evolutions of the state components of the true 
system. We also need to specify the initial conditions for Eq. (15) for the internal variables of the form < ... >^. 
For our experiment, these are all initially set equal to random numbers drawn from a gaussian distribution with zero 
mean and standard deviation equal to 0.1. v is equal to 10^^/2, so that the sliding exponential time average < ... >i, 
is performed over a time ^ 200 that is long compared to the characteristic time Tg for the evolution of a Rossler 
system. We estimate the latter time to be about Ts ~ 6, which is the average measured time interval between two 
consecutive peaks of xi{t). Numerical results are shown in Figs. 1-2. Figure 1 shows the time evolutions of xi(t) and 
x'i{t) (respectively, X2{t) and X2{t), X3{t) and x'^{t)) at the beginning and at the end of the simulation (note that at 
the end of the run, synchronization is attained for all the three state variables). Figure 2 shows the time evolution 
of some of the estimated parameters (in red), when compared to the corresponding true values in system ([3]) (black 
dotted lines). The values of all the M estimated parameters a'yj.(t), b'^j{t), c'^{t) at the end of the run {t ~ 10*) are 
reported in Table 1. It is seen that both the coefficients that have true values equal zero and those that have true 
values different from zero are accurately estimated. 

We have repeated the experiment above many times, and we have observed either success or failure of our strategy 
depending on the random choice of initial conditions for Eqs. (4) and (15). We can explain the failure of our 
strategy, observed in some experiments, in terms of the transient dynamics towards synchronization. In fact, if during 
the identification process, some of the coefficients a[,f.,b[j^ and c[ are not correctly identified, these may eventually 
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FIG. 1: The figure shows the time evolutions of xi{t),X2{t),X3{t) (in black) and of x'i{t),x'2{t),x'j{t) (in red) at the beginning 
and at the end of the simulation, ly = 10~^/2, 7 = 2. 



assume values that make x'(<) too large and x(t) diverge. Since this divergence typically occurs on a time-scale, which 
is faster than that on which our strategy operates (and the coefficients are updated), our strategy fails at correctly 
identify the coefficients. Here, to solve this problem, we propose to replace Eq. (4), by the following equation. 



x' = F'(x') +r(i/(x) -i?(x')), 



(17) 
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FIG. 2; The figure shows the time evolution of some of the estimated parameters for the Rossler system (in red), when compared 
to their true values (black dotted lines), v — 10~^/2, 7 = 2. 



where F'(x') = [/i(x), ^(x), /„(x)]^ and 



/»(x) = { 



a, if /i(x) > a, 
/i(x), if |/i(x)|<a, 
-a, if /j(x) < -a, 



(18) 



where a is a given constant. To test this proposed alternative scheme, we have performed numerical simulations, 
involving integration of Eq. (17), (13), and (15). For example, by setting a value of a = 10* in Eqs. (17-18), our 
adaptive strategy was always observed to be successful in yielding the correct identification of the parameters. 
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TABLE L Coefficients estimated by our strategy for the Rossler and the Lorenz systems at the end of a run, t = 10*. 



Rossler Lorenz 


Rossler Lorenz 




Rossler 


Lorenz 


a'ni 


0.000 


0.000 


•^211 


0.000 


0.000 


„' 

"311 


u.uuu 


n nnn 
u.uuu 


^^112 


0.000 


0.000 


^212 


0.000 


0.000 


„' 

"312 


n nnn 
u.uuu 


1 nn 
1 .uu 


o-'iis 


0.000 


0.000 


•^213 


0.000 


-1.00 


„/ 

"313 


1 nnn 

i.UUU 


n nnn 
u.uuu 


^122 


0.000 


0.000 


<^222 


0.000 


0.000 


r,' 

"322 


n nnn 
u.uuu 


n nnn 
u.uuu 


'^123 


0.000 


0.000 


'^223 


0.000 


0.000 


n' 

"323 


n nnn 

u.uuu 


n nnn 

u.uuu 


'^133 


0.000 


0.000 


'^233 


0.000 


0.000 


0333 


0.000 


0.000 


b'n 


0.000 


-10.0 


b'21 


1.000 


28.0 


b'31 


-0.001 


0.000 


b'i2 


-1.00 


10.0 


b'22 


0.165 


-1.00 


b'32 


-0.001 


0.000 


b'i3 


-1.00 


0.000 


b'23 


0.000 


0.000 


b'33 


-10.0 


-2.67 


ci 


0.000 


0.000 


C'2 


0.000 


0.000 


c'3 


0.202 


0.000 



B. An experiment with the Lorenz system 

We now consider an example in which the unknown system to be identified is the Lorenz system, described by 
X = -F(x), X = {xi,X2, ...,X3)"^, i.e., n — 3. For the Lorenz system, the coefficients aijk{t), bij{t), and Ci(t) in Eq. (1) 
are zero, except for 0213 — —1, 0312 — 1, bn — —10, bi2 — 10, 621 — 28, 622 = —1, 633 = —8/3, for which the system 
is chaotic. 

Again, we assume that x(0) is initialized at a random state on the chaotic attractor, while the system in (j3|) is 
evolved starting from a perturbed initial condition (jl6p . with pi = 7.93, p2 — 9.01, pa — 8.62. In our test, we use our 
adaptive strategy described in this paper, with the aim of identifying all the unknown system parameters, a^fc, 5y , 
Ci. We set 7 = 10, and u = 10^^, so that the sliding exponential time < ... >^ is performed over a time = 100 
that is long compared to the characteristic time T' s — 0.97 for the evolution of a Lorenz system (defined as the time 
at which the autocorrelation function of xi{t) decays at one half of its value a,t t — 0). Table I shows a comparison of 
all the M = 30 parameters estimated at the end of the nm {t = lO"') to their true values. As can be seen, by using 
our strategy, we are able to correctly identify all the M = 30 unknown parameters of the chaotic Lorenz system. 

C. An experiment with time varying system parameters 

We now present results illustrating the use of our adaptive strategy to dynamically estimate the evolutions of the 
true system parameters when these are slowly time varying; where by 'slowly' we mean that they evolve on a time 
scale which is much longer than the averaging time 

In this section, we consider that x = F{t, x), with x — {xi, X2, Xn)^ and F{t, x) = x), /2(i, x), /„(<, x)]-^. 



11 



where 



^) = X! X! O'i0kit)XjXk + ^ K{t)Xj + Ci{t). 
j=l k=j 3 = 1 



(19) 



We expect our identification strategy to work, provided that the averaging time v^^ <^ 9, where 0_is_the time scale 
on which the true system parameter drift takes place (for a similar use of our adaptive strategy, see 
a numerical test (2), we have modified the Rossler equations as follows. 



14|). To perform 



F(x) = 



-a;2 - X3 
xi + 0.165a;2 
0.2 + (.Ti - 10 + sin(cji))x3 



(20) 



i.e., all the parameters are the same as in Sec. III. A (i.e., time invariant), but 633 (t) = 10 + sin(ti;t). System ([^0]) is 
observed to be chaotic in the range 9 < 633 < 11. 

In our test, we use our adaptive strategy described in this paper, Eqs. (4,13,15), with the aim of identifying all the 
unknown system parameters, aijk{t), bij{t), Ci{t). We set 7 = 2, and h' — 0.1, so that Tg < v^^ = 10 ^ 

For the case lo — (lOOT^)"^, Fig. 3 shows the time evolutions of a'^i^^it) when compared to a3i3(i), and b'^^{t) when 
compared to 633 (i). The estimated parameters are seen to accurately reproduce the evolutions of the true ones. Fig. 
4 shows the time averaged identification error — ^ /lo^-A l^'ssl^) ~ b^i{t)\dt (A — 0.3 x 10*) versus loTs for 7 = 2, 
V = 0.1. The point indicated with a full dot corresponds to ojTs = 0.01 (the case shown in Fig. 3); as can be seen, the 
identification error becomes large if ujTs <C 1 is not satisfied. 

IV. NOISE ANALYSIS 

Here we consider the effects of measurement noise on our identification procedure. Namely we replace Eq. ^ by 



x' = F'(x') + r(iJ(x")-H(x')), 



(21) 



where the noisy state x" = (x", , xJJ)-^ is equal to x" = x + r;^, 77 is a scalar gain and ^ = [S,i,£,2, ■■■,£,mV 
is an m-dimensional vector, whose components are normalized noise terms; we choose ^i(t) = zpiei{t), where for 
each i, pi is the standard deviation of the time evolution of Xi for the true system; ei{t) at each time step of our 
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FIG. 3: The figure shows the comparison between the time evolutions of two estimated parameters for system (|20|l (in red), 
when compared to their true values (black dotted lines), uTa = 0.01, v — 0.1, -y — 2. In the example, the true fe33(t) is a 
function of time, i.e., 633(4) = —10 + sin(ajt). 
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FIG. 4: The figure shows the time averaged identification error E} versus loTs for y — 2, v = 0.1. The point indicated with a 
full dot corresponds to ujTa = 0.01 (the case shown in Fig. 3). 

numerical integration are zero-mean independent random numbers of unit variance drawn from a normal distribution; 
z is normalization factor, which is chosen so that 77 1 makes the noise cause the relevant state component to diffuse 
by an amount that, over a time interval T^, is roughly as big as the amplitude of variation of the relevant chaotic 
state component, and is given by z = \J^^ where r is the time step of our integration method and Ts is the period 
of one oscillation of the true system. 
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Measurement noise introduces a rapid small-scale variation to the measurement time series x"(t), thus making the 
meaning of the time derivative dx"(t)/(it that would appear in Eq. ([5]) questionable. Therefore we replace Eq. ([5]) by 



, . . . , If-, 



(22) 



where we obtain x*{t) by low pass filtering the noisy state x"{t), with cutoff frequency T. Then x*i to be inserted in 
(p2| is obtained from the differential equation 



^*^{t) + ^X*{t)^^xnt), 



(23) 



where T is chosen to satisfy Tg >> T >> t. In our numerical experiment, we choose t = 10~^, and T = 4 x 10""^. 

We monitor the robustness of our identification strategy with respect to increasing values of the noise rj. To this 
aim, we introduce the following identification error measure, 



n n n 



= ]^[E EE K-fe w - ^^^^ + E E %^^) - ^-^^ + E 1^^ w - ^^ii- 



(24) 



i=l j=l k=j 



i=l 3 = 1 



i=l 



We have performed numerical tests to investigate how the time averaged identification error Ej — ^ /loi-A Ei{t)dt 
depends on ry (A = 0.3 x 10"^). Our simulations show that Ei remains small through the range tested, and for example 
is less than 1% in the range < ?y < 1, when v ^ A y. 10^^, 7 = 5. 



V. THE CASE OF INCONSISTENT MODEL EQUATIONS 



In this section, we suppose that the dynamics of the true system deviates from polynomial form (i.e., /i(x) cannot 
be written as a given degree polynomial in x) , but we still do our adaptive procedure using a polynomial model as in 
Eq. (I19p . As a first attempt, we continue to assume that the system dynamics can be approximately modeled as a 
degree 2-polynomial as in ([3]) and in so doing we want to test the limits of this approach. Let us consider for example 
that the true system dynamics is described by x = ^"(x), where -F'(x) is given by. 



i^(x) 



-X2 — a;3 + Kpisech(xi) 
xi + 0.165a;2 
0.2 + (xi - 10)x3 



(25) 
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Note that for k = 0, Eq. p5|) is equivalent to the Rossler Eq. ([2]), while for k ^ the system dynamics includes a 
transcendental function, i. e., a function that cannot be expressed as a finite degree polynomial. In what follows we 
implement the same adaptive strategy described in Sec. II and test its effectiveness for increasing values of k. 

In following our adaptive procedure, because the true and model systems are different, the evolving model pa- 
rameters from our system of adaption equations do not settle into constant values. Rather we observe that they 
time-asymptotically fluctuate around some constant value. Thus we take for our identification of the model parame- 
ters the time average of these quantities. 

In order to formulate a measure of the effectiveness of our identification procedure with an inconsistent model, we 
note that for practical purposes, one is often interested in how well a model is able to reproduce the behavior of the 
real system. In particular, a sensible question to ask would be how well the model equations obtained through our 
adaptive strategy forecast the future behavior of the true system, and, in particular how far in the future are such 
forecasts reliable. In what follows we provide a partial answer to such a question. 

We have performed numerical experiments in which we evaluate the error of the obtained model system when it is 
used to forecast the evolution of the true system. We consider two cases, k = 0.1 and k — 0.2. In Fig. [5] we have 
monitored the forecast error pi~^\xi{t) — x'i{t)\ as function of t, when the model and the system are uncoupled and 
evolved from the same initial condition a;i(0) = ^[{O) (which we have taken to be a randomly chosen point from the 
Rossler attractor). The results in Fig. [5] have been averaged over 500 different choices of the initial conditions. From 
Fig. O we observe that for example, if we set our prediction time to be the length of time over which the prediction 
error remains less than 10%, it is about 12 (about 20) in the case of k = 0.1 (k = 0.2). In particular, the latter time 
is more than three times the characteristic time scale of the true system ~ 6. 
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FIG. 5: The figure shows the prediction error pi~^\xi{t) — x'i{t)\ as function of time for the cases of « = 0.1 and k = 0.2. The 
results have been averaged over 500 different choices of the initial conditions. 
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VI. CONCLUSION 

In this paper, we have introduced a new strategy to identify the parameters of an unknown chaotic dynamical 
system. We aim at synchronizing the real unknown system with another system 'in silico', whose parameters are 
adaptively evolved to converge on those of the real one. 

As a first attempt, we have assumed that the differential equations governing the system dynamics are expressible, 
or approximately expressible, in terms of polynomials of an assigned degree. For this case, our strategy relies on 
the assumption that the only necessary information about the true system is the dimensionality of its state vector 
and the order of the polynomials. Under these conditions, we have shown that we are able to extract the whole 
set of parameters of the unknown system from knowledge of the dynamical evolution of its state vector and its first 
derivative. Our procedure relies on the minimization of appropriately defined potentials that are zero when both 
the system state and parameters are correctly identified. Interestingly, our strategy is effective in detecting which 
parameters are/are not zero, and in obtaining correct estimates for those that are not zero. 

We have further considered the effects of measurement noise and we have proposed an alternative scheme that 
works when only knowledge of the dynamical evolution of true system state vector is available, which has been shown 
to be effective even in the presence of relatively high noise. 

We have also analyzed a situation in which the model fitting function basis is slightly inconsistent with the true 
system dynamics, and, for this case, we have evaluated how well the obtained model is able to forecast the future 
behavior of the true system, i.e., how far into the future it is able to forecast its evolution. 

As a further application, we presented the possibility of extending our approach to the case in which the parameters 
to be identified are slowly varying in time (i.e., on a time scale that is slower than The general strategy can 

also be used if one has a system of known form with several unknown parameters. 

This work was supported by the U.S. Office of Naval Research, contract N00014-07-1-0734. 



[1] J. L. Hindmarsh and R. M. Rose, Proc. R. Soc. London B 221 (1984). 

[2] T. Matsumoto, IEEE Transactions on Circuits and Systems CAS 31, 1055 (1984). 

[3] T. Matsumoto, L. O. Chua, and M. Komuro, IEEE Transactions on Circuits and Systems CAS 32, 798 (1985). 

[4] A. Colicn, B. Ravoori, T. E. Murphy, and R. Roy, Phys. Rev. Lett. (2008). 

[5] J. P. Crutchfield and B. S. McNamara, Complex Systems 1 (1987). 



16 

[6] H. Nijmeijer and I. M. Y. Mareels, IEEE Trans. Circuits Syst. I 44, 882 (1997). 

[7] A. Pogromsky and H. Nijmeijer, International Journal of Bifurcation and Chaos 8, 2243 (1998). 

[8] E. O. Paul So and W. P. Dayawansa, Phys. Rev. E 49, 2650 (1994). 

[9] H. D. I. Abarbanel, D. R. Creveling, and J. M. Jeanne, Phys. Rev. E 77, 016208 (2008). 
[10] D. R. Creveling, P. E. Gill, and H. D. I. Abarbanel, Phys. Lett. A 372, 2640 (2008). 
[11] W. Yu, G. Chen, J. Gao, J. Lu, and U. Parlitz, Phys. Rev. E 75, 067201 (2007). 
[12] M. S. Suarez-Castafion, C. Aguilar-Ibafiez, and F. Flores-Ando, Physics Letters A 317, 265 (2003). 
[13] C. Aguilar-Ibaiiez, E. Hernandez-Rubio, and M. S. Suarez-Castafion, Rivista Mexicana de Fisica 53, 436 (2007). 
[14] F. Sorrentino and E. Ott, Phys. Rev. Lett. 100, 114101 (2008). 
[15] L. Pecora and T. Carroll, Phys. Rev. Lett. 80, 2109 (1998). 



